UNCLASSIFIED 


Defense  Technical  Information  Center 
Compilation  Part  Notice 

ADPO 13722 

TITLE:  Orthogonal  Distance  Fitting  of  Parametric  Curves  and  Surfaces 
DISTRIBUTION:  Approved  for  public  release,  distribution  unlimited 


This  paper  is  part  of  the  following  report: 

TITLE:  Algorithms  For  Approximation  IV.  Proceedings  of  the  2001 
International  Symposium 

To  order  the  complete  compilation  report,  use:  ADA412833 

The  component  part  is  provided  here  to  allow  users  access  to  individually  authored  sections 
of  proceedings,  annals,  symposia,  etc.  However,  the  component  should  be  considered  within 
the  context  of  the  overall  compilation  report  and  not  as  a stand-alone  technical  report. 

The  following  component  part  numbers  comprise  the  compilation  report: 

ADPO  13708  thru  ADPO  13761 


UNCLASSIFIED 


Orthogonal  distance  fitting  of  parametric 
curves  and  surfaces 

Sung  Joon  Ahn,  Engelbert  Westkamper,  and  Wolfgang  Rauh 

Fraunhofer  Institute  for  Manufacturing  Engineering  and  Automation  (IPA) 
Nobelstr.  12,  70569  Stuttgart,  Germany 
{sja;  wke;  wor}@ipa. fhg.de 


Abstract 

Fitting  of  parametric  curves  and  surfaces  to  a set  of  given  data  points  is  a relevant 
subject  in  various  fields  of  science  and  engineering.  In  this  paper,  we  review  the  current 
orthogonal  distance  fitting  algorithms  for  parametric  models  in  a well  organized  and  easily 
understandable  manner,  and  present  a new  algorithm.  Each  of  these  algorithms  estimates 
the  model  parameters  minimizing  the  square  sum  of  the  error  distances  between  the  model 
feature  and  the  given  data  points.  The  model  parameters  are  grouped  and  simultaneously 
estimated  in  terms  of  form,  position,  and  rotation  parameters.  The  form  parameters 
determine  the  shape  of  the  model  feature,  and  the  position/rotation  parameters  describe 
the  rigid  body  motion  of  the  model  feature.  The  new  algorithm  is  applicable  to  any  kind 
of  parametric  curve  and  surface.  We  give  fitting  examples  for  circle,  cylinder,  and  helix 
in  space. 

1 Introduction 

The  use  of  parametric  curves  and  surfaces  is  very  common  and  model  fitting  to  a set  of 
given  data  points  is  a relevant  subject  in  various  fields  of  science  and  engineering.  For 
fitting  of  curves  and  surfaces,  orthogonal  distance  fitting  is  of  primary  concern  because 
of  the  applied  error  definition,  namely  the  shortest  distance  from  the  given  point  to  the 
model  feature  [5,  9].  While  there  are  orthogonal  distance  fitting  algorithms  for  explicit  [3], 
and  implicit  models  [2,  7]  in  the  literature,  we  are  considering  in  this  paper  fitting 
algorithms  for  parametric  models  [4,  6,  8,  10,  11]  (Fig.  1). 

The  goal  of  the  orthogonal  distance  fitting  is  the  estimation  of  the  model  parameters 
minimizing  the  performance  index 

al  = (X  — X')TPTP(X  - X')  (1.1) 

or 

<7o  = dTPTPd , (1.2) 

where  XT  = (Xj, . . . , X£)  and  X,T  = (X'1T, . . . , X'J)  are  the  coordinates  vectors  of  the 
m given  points  and  of  the  m corresponding  points  on  the  model  feature,  respectively. 
Moreover,  dT  = (di, . . . , d,„ ) is  the  distances  vector  with  d,  = ||X;  — X-||,  PTP  is  the 
weighting  matrix.  We  are  calling  the  fitting  algorithms  based  on  the  performance  indexes 
(1.1)  and  (1.2)  coordinate-based  algorithm  and  distance-based  algorithm,  respectively. 


122 


Orthogonal  distance  fitting  of  parametric  models 


Fig.  1.  Parametric  features,  and  the  orthogonal  contacting  point  x-  in  frame  xyz  from 
the  given  point  X*  in  frame  XYZ:  (a)  Curve;  (b)  Surface. 

In  this  paper,  the  model  parameters  a are  grouped  and  simultaneously  estimated 
in  three  categories.  First,  the  form  parameters  ag  (e.g.  three  axis  lengths  a , b,  c of  an 
ellipsoid)  describe  the  shape  of  the  standard  model  feature  defined  in  model  coordinate 
system  xyz  (Fig.  1) 

x = x(ag,u)  with  ag  = (ai, . . . , a;)T  . (1.3) 

The  form  parameters  are  invariant  to  the  rigid  body  motion  of  the  model  feature.  The 
second  and  the  third  parameters  groups,  respectively  the  position  parameters  ap  and  the 
rotation  parameters  ar,  describe  the  rigid  body  motion  of  the  model  feature- in  machine 
coordinate  system  XYZ: 

X = R-1x  + X0  or  x = R(X  — XD),  (1.4) 

where  R = RKRVR„,  = (rx  r2  r3)T , R_1  = RT  , 

ap  = X0  — (XOJ  Yo,  Z"0)  , and  ar  — (ui,  y?,  ac) 

A subproblem  of  the  orthogonal  distance  fitting  of  a parametric  model  is  the  finding 
of  the  location  parameters  {ut}™  j,  which  represent  the  nearest  points  {X'}™  x on  the 
model  feature  from  each  given  point  {Xj}^.  The  model  parameters  a and  the  location 
parameters  {u,}("  j will  generally  be  estimated  through  iteration.  By  the  total  method  [6, 
10],  a and  {ui}'£=1  will  be  simultaneously  determined,  while  they  are  to  be  separately 
estimated  by  the  variable-separation  method  [4,  8, 11]  in  a nested  iteration  scheme.  There 
could  be  four  combinations  for  algorithmic  approaches  as  shown  in  Table  1.  One  of  the 
algorithmic  approaches  in  Table  1 results  in  an  obviously  underdetermined  linear  system 
for  iteration,  thus,  it  has  no  practical  application.  We  describe  and  compare  the  realistic 
three  algorithmic  approaches  in  the  following  sections. 
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Algorithmic  approaches 

Distance-based  algor.  Coordinate-based  algor. 

Total  method 

Variable-separation  method 

Underdetermined  system  I (ETH  [6,  10]) 

II  (NPL  [4,  11])  III  (FhG,  this  paper) 

Tab.  1.  Orthogonal  distance  fitting  algorithms  for  parametric  models. 


2 Orthogonal  distance  fitting  algorithm  I (ETH) 

The  ETH  algorithm  [6,  10]  is  based  on  the  performance  index  (1.1),  and  simultaneously 
estimates  the  model  parameters  a and  the  location  parameters  {u,}™=1  for  the  nearest 
points  on  the  model  feature.  We  introduce  the  new  estimation  parameters  vector  b 
containing  a and  {uj}™  j as  follows, 

bT  = (aT,  = (*£,  aj,  a*  uj,  . . . , <) . 

The  parameters  vector  b minimizing  the  performance  index  (1.1)  can  be  determined  by 
the  Gauss-Newton  method 

Ab  = P(X  - X')|* , bk+1  = bk  + aAb , (2.1) 

k 

with  the  Jacobian  matrices  of  each  point  X'  on  the  model  feature,  from  (1.3)  and  (1.4) 
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A disadvantage  of  the  ETH  algorithm  is  that  the  storage  space  and  the  computing  time 
cost  increase  very  rapidly  with  the  number  of  the  data  points,  unless  the  sparse  linear 
system  (2.1)  is  handled  beforehand  by  a sparse  matrix  algorithm. 

3 Orthogonal  distance  fitting  algorithm  II  (NPL) 


The  NPL  algorithm  [4, 11]  is  based  on  the  performance  index  (1.2),  and  separately  estim- 
ates the  model  parameters  a and  the  location  parameters  {ujj-lLj  in  a nested  iteration 
scheme  . . ■>,,  v,, 

mm  mm  ag({xda,u)}£Li)  • 

a i 

The  inner  iteration  determines  the  location  parameters  {u'} £Lj  for  the  minimum  distance 
points  {X'}]Tj  on  the  current  model  feature  from  each  given  point  {Xj}£Li>  and,  the 
outer  iteration  updates  the  model  parameters.  In  this  paper,  in  order  to  implement  the 
parameters  grouping  of  aT  — (aj,aj,a^),  we  have  modified  the  initial  NPL  algorithm. 


3.1  Orthogonal  contacting  point 

For  each  given  point  x*  = R(X,  — XQ)  in  frame  xyz,  we  determine  the  orthogonal  contact- 
ing point  x'  on  the  standard  model  feature  (1.3).  Then,  the  orthogonal  contacting  point 
X'  in  frame  XYZ  to  the  given  point  X,  will  be  obtained  through  a backward  transform- 
ation of  x(  into  XYZ.  We  are  searching  the  location  parameters  u which  minimizes  the 
error  distance  between  the  given  point  x;  and  the  corresponding  point  x on  the  model 
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feature  (1.3) 


D = (xj  - x(ag,  u))T(xi  - x(ag,  u)) . 


(3.1) 


f(*i,x(a„u))  = i(^)  = -(((J._; 


= 0. 


(3.2) 


The  first  order  necessary  condition  for  a minimum  of  (3.1)  as  a function  of  u is 

-x(ag,u))Tx„^  _ 
x(ag,u))Tx„ 

The  condition  (3.2)  means  that  the  error  vector  (xj— x)  and  the  surface  tangent  vectors 
3x/3u  at  x should  be  orthogonal.  We  solve  (3.2)  for  u by  using  the  Newton  method 
(how  to  derive  the  Jacobian  matrix  df/du  is  shown  in  Section  4). 

df 

— Au  = -f(u)|fc,  Ufc+1  = ufc  + aAu . (3.3) 

3.2  Orthogonal  distance  fitting 

We  update  the  model  parameters  a minimizing  the  performance  index  (1.2)  by  using 
the  Gauss-Newton  method  (outer  iteration) 

dd 

P—  Aa=— Pd|fc,  afc+i  = afe  + aAa . 

da  k 

Prom  dj  = ||Xj  - X'||,  and  equations  (1.3)  and  (1.4),  we  derive  the  Jacobian  matrices  of 
each  orthogonal  distance  di 

ddi  (Xi  - X')T  8X I 


|X 

,T 


XIII  da 


- (X»-XD  Ap-i  /dx  3xdu\  8R  1 
“ HXi-X'll  V {da^dudaj^  da 


dXo 

da 


With  (1.4)  and  (3.2)  at  u — u(, 
(X. 


xifR-S 


>\T 


and  Jdi.. 


(X,  - X 

' HXi-XIl 


R 


du 
! 9x 


= (Xi 


,/\T 


da„ 


<9x 

du 

dR- 

dar 


= (P 


-x. 


is  the  resultant  Jacobian  matrix  for  di.  A drawback  of  the  NPL  algorithm  is  that  the 
convergence  and  the  accuracy  of  3D-curve  fitting  (e.g.  fitting  of  a circle  in  space)  are 
relatively  poor.  2D-curve  fitting  or  surface  fitting  with  the  NPL  algorithm  do  not  suffer 
from  such  problems. 

4 Orthogonal  distance  fitting  algorithm  III  (FhG) 

At  the  Fraunhofer  Institute  IPA  (FhG-IPA),  a new  orthogonal  distance  fitting  algorithm 
for  parametric  models  is  developed,  which  minimizes  the  performance  index  (1.1)  in  a 
nested  iteration  scheme  (variable-separation  method).  The  new  algorithm  is  a general- 
ized extension  of  an  orthogonal  distance  fitting  algorithm  for  implicit  plane  curves  [1]. 
Interested  readers  are  referred  to  [2]  for  the  orthogonal  distance  fitting  of  implicit  sur- 
faces and  plane  curves.  The  location  parameter  values  (u Il-Tj  for  the  minimum  distance 
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points  on  the  current  model  feature  from  each  given  point  {X,}’^,  are  to  be 

found  by  the  algorithm  described  in  Section  3.1  (inner  iteration).  In  this  section,  we 
intend  to  describe  the  outer  iteration  which  updates  the  model  parameters  a minimizing 
the  performance  index  (1.1)  by  using  the  Gauss-Newton  method 


,5X; 

da 


Aa  = P(X  - X') 


ak+l  = a*.  + aAa , 


(4.1) 


with  the  Jacobian  matrices  of  each  orthogonal  distance  point  X',  from  (1.3)  and  (1.4) 
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(4.2) 


The  derivative  matrix  5u/5a  at  u = u'  in  (4.2)  describes  the  variational  behavior  of 
the  location  parameters  u'  for  the  orthogonal  contacting  point  x'  in  frame  xyz  relative 
to  the  differential  changes  of  the  parameters  vector  a.  Purposefully,  we  derive  5u/5a 
from  the  condition  (3.2).  Because  (3.2)  has  an  implicit  form,  its  derivatives  lead  to 


5f  5u  5f  dx(  5f  _ 5f  5u  _ / 5f  5x,  5f\ 

5u  5a  5xj  5a  5a  °r  5u  5a  \ 5xj  5a  ^ 5a  / ’ 

where  5xj/5a  is,  from  x,  = R(X,  - XQ), 


5X; 
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^(X,-X0)-R 


5Xo 
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-R 


— fX-X 
5 (X,  X0 


The  other  three  matrices  5f/5u,  5f/5xi,  and  5f/5a  in  (3.3)  and  (4.3)  are  to  be  directly 
derived  from  (3.2).  The  elements  of  these  three  matrices  are  composed  of  simple  linear 
combinations  of  components  of  the  error  vector  (x,  - x)  with  elements  of  the  following 
three  vector/matrices  5x/5u,  H,  and  G (XHG  matrix): 


Now  (4.3)  can  be  solved  for  5u/5a  at  u = u',  and  the  Jacobian  matrix  (4.2)  and  the 
linear  system  (4.1)  can  be  completed  and  solved  for  the  parameter  update  Aa. 

We  would  like  to  stress  that  only  the  standard  model  equation  (1.3),  without  involve- 
ment of  the  position/rotation  parameters,  is  required  in  (4.4).  The  overall  structure  of 
the  FhG  algorithm  remains  unchanged  for  all  dimensional  fitting  problems  of  parametric 
models.  All  that  is  necessary  for  a new  parametric  model  is  to  derive  the  XHG  matrix 
of  (4.4)  from  (1.3)  of  the  new  model  feature,  and  to  supply  a proper  set  of  initial  para- 
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Fig.  2.  Information  flow  with  the  FhG  algorithm. 


X 5 655320  -1  -1  0 3 4 7 9 

Y 1 3 4 6 5 4 2 0 -2  -5  -7  -8  -10  -9 

Z -3  -1  1 3 5 7 9 11  11  11  11  11  11  10 

Tab.  2.  Fourteen  coordinate  triples  representing  a helix. 

meter  values  a0  for  iteration  (4.1).  An  overall  schematic  information  flow  with  the  FhG 
algorithm  is  shown  in  Fig.  2.  The  FhG  algorithm  shows  robust  and  fast  convergence 
with  2D/3D-curve  and  surface  fitting.  The  storage  space  and  computing  time  cost  are 
proportional  to  the  number  of  data  points.  A disadvantage  of  the  FhG  algorithm  is  that 
it  additionally  requires  the  second  derivatives  <92x/<9ag<9u  as  shown  in  (4.4). 

As  a fitting  example,  we  show  the  orthogonal  distance  fitting  of  a helix.  The  standard 
model  feature  (1.3)  of  a helix  in  frame  xyz  can  be  described  as  follows.  x(ag,  u)  = 
x(r,h,u)  = (rcosu,  rsinu,  hu/2ir)T , with  a constraint  on  the  position  and  rotation 

parameters  

/c(ap,ar)  = (X0  - X)Tr3(w,<£>)  = 0, 

where  r and  h are  respectively  the  radius  and  elevation  of  a helix.  X is  the  gravitational 
center  of  the  given  points  set  and  r3  (see  (1.4))  is  the  vector  of  direction  cosines  of 
the  z-axis.  We  have  obtained  the  initial  parameter  values  from  a 3D-circle  fitting,  and 
a cylinder  fitting,  successively.  The  helix  fitting  to  the  points  set  in  Table  2 with  the 
initial  values  of  h — 10  and  k = -k  terminated  after  0.22s,  8 iteration  cycles  for  ||Aa||  = 
3.2  xlO-7  with  a Pentium  133  MHz  PC  (Table  3,  Fig.  3).  They  were  0.33s,  10  iteration 
cycles  for  ||Aa||=3.6xlO~7  with  the  ETH  algorithm,  and,  1.05s,  61  iteration  cycles  for 
1| Aa||  =8.8xl0~7  with  the  NPL  algorithm.  The  computing  cost  with  the  ETH  algorithm 
increases  rapidly  with  the  number  of  the  data  points.  The  NPL  algorithm  showed  slow 
convergences  with  the  3D-circle  and  the  helix  fitting  (3D-curve  fitting). 
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Parameters  a 

<A) 

r 

h 

x0 

3D-Circle 

5.8913 

8.3850 

— 

5.6999 

a (a) 

— 

0.7355 

— 

0.9939 

Cylinder 

1.6925 

8.2835 

— 

4.7596 

<r(a) 

— 

0.2738 

— 

0.7465 

Helix 

2.2301 

6.1368 

19.5811 

3.8909 

<r(a) 

— 

0.4238 

1.3214 

0.5488 

Zo 

U) 

H> 

K 

-2.7923 

5.2333 

-0.6833 

0.7882 



0.8421 

0.8821 

0.1177 

0.1375 

— 

-3.0042 

4.5081 

-0.4576 

1.1327 

— 

0.4525 

0.6513 

0.3049 

0.2116 

— 

-1.5560 

6.4871 

0.3003 

0.5114 

2.4602 

0.3934 

0.7500 

0.0880 

0.0663 

0.2881 

Tab.  3.  Results  of  the  orthogonal  distance  fitting  to  the  points  set  in  Table  2. 


b Iteration  number 

(a)  (b) 

Fig.  3.  Orthogonal  distance  fitting  to  the  points  set  in  Table  2:  (a)  Helix  fit;  (b)  Con- 
vergence of  the  fit.  Iteration  number  0-3:  3D-circle,  4-12:  circular  cylinder,  and  13-: 
helix  fit  with  the  initial  value  of  h = 10  and  n = n. 

5 Summary 

In  this  paper,  we  have  reviewed  the  current  orthogonal  distance  fitting  algorithms  for 
parametric  curves  and  surfaces  in  an  easily  understandable  manner,  and  presented  a new 
algorithm.  By  each  of  the  algorithms  the  model  parameters  are  grouped  and  simultan- 
eously estimated  in  terms  of  form/position/rotation  parameters.  The  ETH  algorithm  de- 
mands a large  amount  of  storage  space  and  high  computing  cost,  and  the  NPL  algorithm 
shows  relatively  poor  performance  with  3D-curve  fitting.  The  new  algorithm,  the  FhG 
algorithm,  has  no  such  drawbacks  of  the  ETH  algorithm  or  of  the  NPL  algorithm.  A 
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disadvantage  of  the  FhG  algorithm  is  that  it  requires  the  second  derivatives  d2x./dagdu. 
The  FhG  algorithm  does  not  require  a necessarily  good  set  of  initial  parameter  val- 
ues, which  could  also  be  internally  supplied  as  demonstrated  with  the  fitting  examples. 
From  the  viewpoint  of  implementation  and  application  to  a new  model  feature,  the  FhG 
algorithm  is  universal  and  very  efficient.  Merely  the  standard  model  equation  (1.3)  of 
the  new  model  feature  is  eventually  required,  which  has  only  few  form  parameters.  The 
functional  interpretation  and  treatment  of  the  position/rotation  parameters  are  basic- 
ally identical  for  all  parametric  models.  The  storage  space  and  the  computing  time  cost 
are  proportional  to  the  number  of  given  data  points.  Together  with  other  orthogonal 
distance  fitting  algorithms  for  implicit  models  [2],  the  FhG  algorithm  is  certified  by  the 
German  federal  authority  PTB  [5,  9],  with  a certification  grade  that  the  parameter  es- 
timation accuracy  is  higher  than  0.1pm  for  length  unit,  and  O.lprad  for  angle  unit  for 
all  parameters  of  all  tested  model  features. 
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